package fem2.element;

import inf.jlinalg.ArrayVector;
import inf.jlinalg.IMatrix;
import inf.jlinalg.IMatrixReadOnly;
import inf.jlinalg.IVector;
import fem2.MaterialModel;
import fem2.MeshPart;

public class NonlinearStructuralElement extends StructuralElement {

	public NonlinearStructuralElement(MeshPart mp, MaterialModel mm) {
		super(mp, mm);
	}

	/**
	 * Constructs B matrix for nonlinear by summation B-linear and nonlinear
	 * w.r.t physical coordinates
	 * 
	 * <pre>
	 * | (dN1/dX1)+u11*(dN1/dX1)                                u21*(dN1/dX1)                 ... | 
	 * |     u21*(dN1/dX2)                                  (dN1/dX2)+u22*(dN1/dX2)           ... | 
	 * | (dN1/dX2)+u12*(dN1/dX1)+u11*(dN1/dX2)     (dN1/dX1)+u22*(dN1/dX1)+u21*(dN1/dX2)      ... |
	 * 
	 * </pre>
	 * 
	 * 
	 * <blockquote> reference: Meschke, Advanced FEM lecture note part 1
	 * </blockquote>
	 * 
	 * @param xi
	 *            natural coordinate
	 * @return matrix of gradients
	 */
	@Override
	public void computeB(MeshPart mp, IMatrix b, double... xi) {
		double[][] gGradients = mp.gGradientsAt(xi);
		int n = mp.getBasis().getDimension();
		double[][] ju = mp.gJacUAt(xi);
		int dim = mp.getDimension();
		if (dim == 2) {
			double u11 = ju[0][0];
			double u12 = ju[1][0];
			double u21 = ju[0][1];
			double u22 = ju[1][1];

			for (int i = 0; i < n; i++) {
				double gra1 = gGradients[i][0];
				double gra2 = gGradients[i][1];

				b.set(0, 2 * i, gra1 + u11 * gra1);
				b.set(0, 2 * i + 1, u21 * gra1);
				b.set(1, 2 * i, u12 * gra2);
				b.set(1, 2 * i + 1, gra2 + u22 * gra2);
				b.set(2, 2 * i, gra2 + u12 * gra1 + u11 * gra2);
				b.set(2, 2 * i + 1, gra1 + u22 * gra1 + u21 * gra2);
			}
		} else if (dim == 3) {
			double u11 = ju[0][0];
			double u12 = ju[1][0];
			double u13 = ju[2][0];
			double u21 = ju[0][1];
			double u22 = ju[1][1];
			double u23 = ju[2][1];
			double u31 = ju[0][2];
			double u32 = ju[1][2];
			double u33 = ju[2][2];

			for (int i = 0; i < n; i++) {
				double gra1 = gGradients[i][0];
				double gra2 = gGradients[i][1];
				double gra3 = gGradients[i][2];

				b.set(0, 3 * i, gra1 + u11 * gra1);
				b.set(0, 3 * i + 1, u21 * gra1);
				b.set(0, 3 * i + 2, u31 * gra1);

				b.set(1, 3 * i, u12 * gra2);
				b.set(1, 3 * i + 1, gra2 + u22 * gra2);
				b.set(1, 3 * i + 2, u32 * gra2);

				b.set(2, 3 * i, u13 * gra3);
				b.set(2, 3 * i + 1, u23 * gra3);
				b.set(2, 3 * i + 2, gra3 + u33 * gra3);

				b.set(3, 3 * i, gra2 + u12 * gra1 + u11 * gra2);
				b.set(3, 3 * i + 1, gra1 + u22 * gra1 + u21 * gra2);
				b.set(3, 3 * i + 2, u32 * gra1 + u31 * gra2);

				b.set(4, 3 * i, u13 * gra2 + u12 * gra3);
				b.set(4, 3 * i + 1, gra3 + u23 * gra2 + u22 * gra3);
				b.set(4, 3 * i + 2, gra2 + u33 * gra2 + u32 * gra3);

				b.set(5, 3 * i, gra3 + u13 * gra1 + u11 * gra3);
				b.set(5, 3 * i + 1, u23 * gra1 + u21 * gra3);
				b.set(5, 3 * i + 2, gra1 + u33 * gra1 + u31 * gra3);
			}
		}
	}

	@Override
	public IMatrixReadOnly computeKte(IMatrix kte) {
		kte = (IMatrix) computeKge(kte);
		kte = (IMatrix) computeMaterialKe(kte);
		return kte;
	}

	/**
	 * compute material stiffness matrix
	 * 
	 * @param ke
	 * @return
	 */
	public IMatrixReadOnly computeMaterialKe(IMatrix ke) {
		super.computeKte(ke);
		return ke;
	}

	/**
	 * <blockquote>refer to equation (1.32)~(1.33), Meschke's advanced FEM
	 * lecture note </blockquote>
	 */
	@Override
	public IVector computeRegularStrain(double... xi) {
		MeshPart mp = getMeshPart();
		double[][] ju = mp.gJacUAt(xi);
		int dim = mp.getDimension();
		int ntrs = dim * (dim + 1) / 2;
		IVector e = new ArrayVector(ntrs);
		if (dim == 2) {
			double u11 = ju[0][0];
			double u12 = ju[1][0];
			double u21 = ju[0][1];
			double u22 = ju[1][1];
			e.set(0, 0.5 * (2 * u11 + u11 * u11 + u21 * u21));
			e.set(1, 0.5 * (2 * u22 + u22 * u22 + u12 * u12));
			e.set(2, u12 + u11 * u12 + u21 + u21 * u22);
		} else if (dim == 3) {
			double u11 = ju[0][0];
			double u12 = ju[1][0];
			double u13 = ju[2][0];
			double u21 = ju[0][1];
			double u22 = ju[1][1];
			double u23 = ju[2][1];
			double u31 = ju[0][2];
			double u32 = ju[1][2];
			double u33 = ju[2][2];
			e.set(0, 0.5 * (2 * u11 + u11 * u11 + u21 * u21 + u31 * u31));
			e.set(1, 0.5 * (2 * u22 + u12 * u12 + u22 * u22 + u32 * u32));
			e.set(2, 0.5 * (2 * u33 + u13 * u13 + u23 * u23 + u33 * u33));
			e.set(3, u12 + u21 + u11 * u12 + u21 * u22 + u31 * u32);
			e.set(4, u23 + u32 + u12 * u13 + u22 * u23 + u32 * u33);
			e.set(5, u13 + u31 + u11 * u13 + u21 * u23 + u31 * u33);
		}
		return e;
	}

}
